Changes in Behavior and Neural Dynamics across Adolescent Development

Adolescence is an important developmental period, during which substantial changes occur in brain function and behavior. Several aspects of executive function, including response inhibition, improve during this period. Correspondingly, structural imaging studies have documented consistent decreases in cortical and subcortical gray matter volume, and postmortem histologic studies have found substantial (∼40%) decreases in excitatory synapses in prefrontal cortex. Recent computational modeling work suggests that the change in synaptic density underlie improvements in task performance. These models also predict changes in neural dynamics related to the depth of attractor basins, where deeper basins can underlie better task performance. In this study, we analyzed task-related neural dynamics in a large cohort of longitudinally followed subjects (male and female) spanning early to late adolescence. We found that age correlated positively with behavioral performance in the Eriksen Flanker task. Older subjects were also characterized by deeper attractor basins around task related evoked EEG potentials during specific cognitive operations. Thus, consistent with computational models examining the effects of excitatory synaptic pruning, older adolescents showed stronger attractor dynamics during task performance. SIGNIFICANCE STATEMENT There are well-documented changes in brain and behavior during adolescent development. However, there are few mechanistic theories that link changes in the brain to changes in behavior. Here, we tested a hypothesis, put forward on the basis of computational modeling, that pruning of excitatory synapses in cortex during adolescence changes neural dynamics. We found, consistent with the hypothesis, that variability around event-related potentials shows faster decay dynamics in older adolescent subjects. The faster decay dynamics are consistent with the hypothesis that synaptic pruning during adolescent development leads to stronger attractor basins in task-related neural activity.

Histologic analysis of postmortem tissue has also shown that excitatory synapses in cortex, which primarily occur on spines, peak between ages 5 and 10, and then decrease exponentially into the mid-20s (Huttenlocher, 1979;Huttenlocher and Dabholkar, 1997;Petanjek et al., 2011).Some studies have found that up to 40% of synapses are pruned during adolescence (Petanjek et al., 2011).Furthermore, synaptic density decreases earlier in primary sensory areas than in association areas (Huttenlocher and Dabholkar, 1997;Elston et al., 2009).Examination of resting-state oscillations in EEG and MEG have also found that low-frequency oscillations decrease during adolescence, and high-frequency oscillations increase (Whitford et al., 2007;Segalowitz et al., 2010;Candelaria-Cook et al., 2022;Rempe et al., 2023).These changes in oscillations correlate with developmental changes in gray matter and white matter (Candelaria-Cook et al., 2022).Evoked responses also develop over adolescence, with early sensory-motor peaks decreasing in amplitude and in peak latency.Other evoked responses associated with cognitive control, performance monitoring and working memory, change in both morphology and location, like the P3 response seen in go/no-go tasks and the error-related negativity (ERN) response, both increasing in amplitude over the adolescent years (Jonkman, 2006;Segalowitz et al., 2010;Tamnes et al., 2013).However, few studies use measures of brain function to test mechanistic theories related to synaptic pruning, and the current report provides such mechanistic data.
Behavior undergoes substantial change during adolescence in ways that could relate to pruning.Although basic perceptual abilities reach adult levels by late childhood (Leat et al., 2009), executive processing, such as working memory and inhibitory control, depend on prefrontal cortex (Goldman- Rakic, 1996) and continue to develop into adulthood (Larsen and Luna, 2018;Theodoraki et al., 2020;Ferguson et al., 2021).Reinforcement learning (RL), which is the ability to learn to select actions that are rewarding and avoid actions that are harmful (B.B.Averbeck and Costa, 2017;B. Averbeck and O'Doherty, 2022), also improves (Lange-Küttner et al., 2012, 2021;Nussenbaum and Hartley, 2019).While RL is often thought to be driven by subcortical plasticity, recent work suggests that working memory may also play an important role in RL (Collins and Frank, 2012;Yoo and Collins, 2022).Cognitive control, measured with speeded response tasks like the Eriksen Flanker task, also improves during adolescence (Luna et al., 2010).
Few mechanistic theories link changes in behavior during adolescence with changes in the brain.It has often been assumed that maturation of prefrontal cortex influences behavior (Caballero et al., 2016).This, however, leaves open the question of mechanisms that link prefrontal cortex maturation to changes in behavior.Recently, a computational model was developed that explicitly linked pruning in recurrent cortical networks to improvements in behavior (B.B.Averbeck, 2022).This study trained recurrent neural networks on working memory and RL tasks.It was found that, if recurrent connections in the network were pruned during training, networks performed better than networks that were given the same amount of training, but not pruned.The current study tests this theory using longitudinal data in children and adolescents.
In simulated data, pruning led to working memory networks that were more resistant to distraction, and RL networks that learned to select rewarding options more consistently, similar to what is seen during adolescent development (Nussenbaum and Hartley, 2019).Furthermore, the pruning led to a characteristic signature in the neural dynamics.Specifically, a metric known as the Lyapunov exponent decreased with pruning.These decreases relate to changes in attractor dynamics and resistance to network perturbation.Thus, pruned neural networks have smaller Lyapunov exponents and stronger attractor dynamics than unpruned networks.Although the network models were trained to reproduce behavior, the pruning led to consistent changes in the neural dynamics, and it is the changes in the neural dynamics that we examine in the present manuscript.
We sought to test the hypothesis that adolescent development, and the corresponding pruning that takes place, leads to stronger attractors during task performance, and therefore measurable changes in dynamics.To do so, we fit dynamical systems models to EEG data collected longitudinally in a large cohort of human subjects, at 12, 15, and 18 years of age.At each age point the subjects performed the Eriksen Flanker task while EEG data were collected.Consistent with model predictions, we found that time-varying linear estimates of the neural dynamics were consistent with stronger dynamical attractors in older subjects.We also found that these changes in dynamics correlated with improvements in reaction time.

Data
Children were recruited for a larger longitudinal study at the University of Maryland assessing behavior (Fox et al., 2001) and EEG markers of error response over development (Buzzell et al., 2017;Filippi et al., 2020).Participants included in this analysis were selected as they completed a Flanker task with EEG at age 12, 15, and/or 18 years.Not all subjects were sampled at all time points.
We preprocessed a total of 484 EEG recordings from 213 participants.We included time locked data with at least 50 trials passing our quality assessment criteria giving us a clean sample of 324 stimulus locked and 307 cue locked and response locked datasets.The clean time locked datasets came from 331 separate recording sessions from 179 children (104 female, 75 male).The average ages at the three recording time points were 13.1 (SD 0.5), 16.0 (SD 0.5), and 18.4 (SD 0.6) for the 12-, 15-, and 18-year target age groups, respectively.
The task consisted of 384 trials separated into 12 blocks.At the end of each block, the computer presented feedback based on the child's performance.In blocks where accuracy was 90% or above, the feedback was "Respond faster."In blocks where accuracy was between 75% and 90%, the feedback was "Good job."And in blocks where accuracy was below 75%, the feedback was "Be more accurate."This feedback was designed so that children produced a sufficient number of errors.
Each trial consisted of presentation of a fixation cross visual cue (lasting ;300-600 ms), followed by a visual stimulus composed of five horizontally aligned arrowheads presented for 200 ms.The stimulus was followed by a blank screen for ;1900 ms.
Participants were instructed to indicate as quickly as possible with a button press the direction of the central arrow.The task had four stimulus types appearing with equal probability: two congruent stimuli where all five arrows pointed in the same direction (n(, o)), and two incongruent stimuli with the "flanking" arrows pointing in the opposite direction from the central arrow ((.(, ),)).
A 128-channel HydroCel Geodesic Sensor Net and EGI software (Electrical Geodesic) were used for all data collection.Data were sampled online at 250 Hz and referenced to the vertex electrode Cz.Data were converted to EEGLAB format before preprocessing trough the MADE pipeline (Debnath et al., 2020).Preprocessing included deleting the outer cap electrodes with the largest artefacts (going from 128 to 104 electrodes), eliminating 1s data segments with noise spikes, high-pass filtering at 0.2 Hz, lowpass filtering at 50 Hz, cleaning of eye motion artefacts with Independent Components Analysis (ICA) with the ADJUST toolbox (Mognon et al., 2011), segmenting the data with reference to the cue, stimulus, and motor response, and finally referencing all sensors from Cz to the electrode average.The full preprocessing code is available on GitHub, and all preprocessing and subsequent analyses were conducted with MATLAB 2020a (The MathWorks) and EEGLAB toolbox v.2022.0(Delorme and Makeig, 2004).

Analysis
The event-related potential (ERP) spectral density is dominated by low frequencies (,10 Hz).To focus on the evolution of the evoked response components we low pass filtered the EEG data at 15 Hz and down sampled to 30 Hz.All trials in which participants responded to the stimulus were included (i.e., all correct and commission error trials) in both congruent and incongruent conditions.
Our goal was to fit linear dynamical systems models to examine the dynamical properties of the EEG signals around the mean responses.We considered four task conditions: congruent versus incongruent, and correct versus commission errors, crossed.There was no evidence that left/ right responses were coded in the EEG signal (p .0.05) so we did not break-out the response type as separate conditions.To reduce the spatial dimensionality of our EEG data we applied principal components analysis (PCA).We calculated the PCA coefficients on the time concatenated stimulus locked and response locked conditions (on the grand-averaged data from all age groups).The PCs were estimated by computing a covariance matrix across these vectors such that the dimensionality of the covariance matrix was 104 Â 104 (given by the number of sensors).We then extracted PCs from this matrix.We selected the first two components for the analysis of the dynamics.By extracting the PCs from the average ERPs, we were targeting a low-dimensional representation in the space of the task-related signals.Note that the PCs have no a-priori theoretical relationship with our task, or with the predictions of the computational model.They were used as a signal processing tool to identify a low-dimensional space of task-related evoked activity in which we could test our hypothesis.
In the next step we projected the data from each time point in each trial into the PC space.This reduced the dimensionality from the 104 sensors to two PC dimensions.We then z-scored each PC for each subject (Fig. 1A).The z-score was computed by first computing the mean and SD across all time points for each component, and then using these to standardize the data.Thus, we subtracted one overall mean from each channel and divided by the overall SD.The resultant time-series in each PC had a mean of 0 and a SD of 1 (Fig. 1B).
Data were next segmented in windows around the presentation of the fixation cue (À600-500 ms relative to the cue), presentation of the flanker stimuli (À400-600 ms relative to the flanker stimuli), and the subject's response (À500-500 ms relative to button responses) for each trial.The dynamics analysis characterizes the temporal evolution (i.e., relaxation) of fluctuations around the mean evoked activity sðtÞ, where t is time in trial.We assessed the relaxation back to the mean response sðtÞ by estimating the decay of the single trial residuals in the 2D PCA space.The residual (Fig. 1C,D), x i;j ðtÞ, is equal to the difference between the single trial activity s i;j ðtÞ, where i is trial and j is condition, and the average across trials for that condition s j ðtÞ, i.e., x i;j ðtÞ ¼ s i;j ðtÞ À s j ðtÞ.Our analyses were then fit to this residual activity (Galgali et al., 2023).The dynamical systems models characterize the rate at which noise in the residual activity decays back to the mean.
We next fit an autoregressive (AR1), linear model to the 2D residual neural activity, x t ð Þ, where we have dropped the subscripts for trial and condition to simplify notation: The variable « ðtÞ is the noise term, assumed to be independent identically distributed (i.i.d.) Gaussian, and A is the 2 Â 2 dynamics matrix, where the two dimensions correspond to the two PCs.The variable x t ð Þ is correspondingly a 2D vector of the residual activity at time t.The model was fit to data points in a 100 ms moving window (at 30 Hz, this corresponds to three time points from each trial sampled every ;33 ms).Thus, data points in a 100-ms window, relative to the trial event being considered (e.g., the cue), from all trials were pooled and used to fit the model.This resulted in a time-dependent estimate of the matrix A. We then extracted the eigenvalues of matrix A and subjected them to further random effects analyses to assess our hypothesis.Eigenvalues of A closer to 0 indicate a faster decay and eigenvalues near 1 would correspond to a slower decay.We expect the eigenvalues of A to be between 0 and 1, and to be smaller in older adolescents.

Statistical tests
According to our hypotheses we expect the eigenvalues of the matrix A to show modulation with the behavioral task and to have smaller values in older adolescents (i.e., a negative relationship between age and eigenvalues).We first tested the eigenvalues with the analysis window locked to each task condition (cue, stimulus, and response) for effects of age.To account for multiple EEG recordings belonging to the same children at different ages we used a linear mixed effect model with a random effect of participant: We conducted a similar analysis to examine the effect of sex with and without age: We also tested the hypothesis that the eigenvalues of A would relate to task performance, with smaller eigenvalues being related to shorter reaction times (positive relationship) and higher task accuracy (negative relationship).Again, to account for the longitudinal recordings we tested these hypotheses with linear mixed effects models: Mixed effects models were run independently for each eigenvalue and time window.To account for multiple comparisons over time, we applied a one dimensional temporal clustering model with the threshold free cluster enhancement algorithm on the fixed effects t statistics (Smith and Nichols, 2009).The null distribution for the temporal cluster was obtained by bootstrap sampling the fixed effect variable for all time points, obtaining the t statistics of the mixed effect model on the bootstrapped time course, applying temporal clustering and picking the highest (or lowest, depending on test hypothesis) cluster value for that time course.The bootstrap procedure was repeated 10,000 times for each hypothesis test to obtain a null distribution for the temporal clusters.For the tests of age effects on eigenvalues we conducted tests locked to three different task epochs (cue, stimulus, and response) and 2 eigenvalues and specific hypotheses on the direction of the relationships.Therefore, we selected one-tailed Bonferroni corrected significance tests with a-level of 0.05/ 6 ¼ 0.0083.For the analysis of behavioral effects, we conducted nine tests on 2 eigenvalues.Therefore, we selected one-tailed Bonferroni corrected significance tests with a-level of 0.05/18 ¼ 0.0027.

Analysis of time-jittered datasets
To generate datasets with jittered ERPs in single trials, we generated single session data, where each trial's evoked response s b i;j t ð Þ, where the superscript indicates jittered data, was given by: For each trial, t was sampled from a mean 0 Gaussian distribution where the SD of the Gaussian was given by the average SD for the corresponding age group.The variable « ðtÞ indicates zero mean Gaussian noise with a SD given by the average across all subjects and groups.

Results
We conducted 128 channel EEG recordings in adolescents while they executed the flanker task.Data were collected at three ages longitudinally for a total of 331 recordings from 179 adolescents (after quality control).In the task, subjects executed 384 trials in 12 blocks in a single session.Subjects were first shown a cue (plus sign), and after a variable delay, they were given a stimulus (Fig. 2A).Their task was to indicate with a key press whether the central arrow pointed left or right.The task had congruent and incongruent conditions.In the congruent condition all arrows pointed in the same direction, and in the incongruent condition the flanking arrows pointed in the opposite direction of the central arrow.At the end of each block, subjects were given feedback which was "good job," "be more accurate," or "go faster," to keep their accuracy between 75% and 90% correct.We found that participants were more accurate [congruent accuracy 0.97 (SD 0.04), incongruent 0.77 (SD 0.10), paired-sample t (330) ¼ 41.3, p , 0.001], and faster [RT congruent 354 ms (SD 40 ms), incongruent 404 ms (SD 52 ms), t (330) ¼ 42.2, p , 0.001] in congruent than incongruent trials.When we examined the effect of age on behavior, we found that accuracy increased with age [Fig.2B; F (1,329) ¼ 54.72, p , 0.001), and reaction times decreased with age (Fig. 2C; F (1,329) ¼ 35.32, p , 0.001].Thus, subjects became more proficient at the task over adolescent development.We also found a small effect of sex on task performance with males responding faster [males RT 369 ms (SD 41 ms), females 387 ms (SD 47 ms), unpaired-sample t (329) ¼ 3.6, p , 0.001)] and with slightly less accuracy [males 0.86 (SD 0.06), females 0.88 (SD 0.05), unpaired-sample t (329) ¼ 2.5, p ¼ 0.01)].
We began the analysis of the EEG data by calculating a set of principal components in sensor space.The components were based on the average response in each condition (concatenated), with the average taken across all trials and subjects (Fig. 3; see Materials and Methods).The first component was a frontal and occipital component explaining 79% of the variance, while the second was a central/parietal and temporal component explaining 14% of the variance (Fig. 3A).The second component also captured most of the differences between correct and error responses.Component 1 reflected an anterior-posterior gradient and component 2 was more centrally focused (Fig. 3B).The first PC shows a rapid negative response, at around 100 ms, and may reflect visual inputs (Fig. 3C).The second PC separates correct responses from commission errors (i.e., an incorrect button press), especially a few hundred milliseconds after the response (Fig. 3C).
To begin the analysis of the dynamics, the data for individual trials in which subjects made a response were projected into the first two PCs (see Materials and Methods; Fig. 1).The data in each PC was then z-scored within each subject using the mean and SD across all trials and time points.This standardized the variance across subjects and age groups.We then computed the time-varying mean of the z-transformed data for each condition and subtracted the mean for each condition for all trials of that condition.This resulted in the residual activity, around the mean, for each single trial.This was done for both PCs.
We then fit a 2D linear dynamical system model to the residual neural activity in the two PC dimensions in 100-ms moving windows (Materials and Methods).The eigenvalues of the linear dynamical system, extracted from the matrix A (Materials and Methods) characterize the dynamics of the neural activity, around the mean neural activity.The first eigenvalue primarily reflected the first PC and the second eigenvalue primarily reflected the second PC, except around 375 ms after stimulus onset.Small eigenvalues indicate fast dynamics, in which variability in neural activity rapidly relaxes back to the mean, and larger eigenvalues indicate slower dynamics, in which variability in neural activity relaxes slowly back to the mean.A value of 1 would be a perfect integrator, a condition in which the activity would integrate variability instead of relaxing back to the mean, as in a diffusion process.Our hypothesis is that adolescent development will lead to faster dynamics, and therefore eigenvalues will get smaller with age.This follows from computational modeling in recurrent neural networks, looking at the effects of pruning on working memory and reinforcement learning (B.B.Averbeck, 2022).
To test this hypothesis, we examined age-dependent changes in the time-dependent eigenvalues.The dynamics models were fit in 100-ms windows time locked to cue onset, stimulus onset, or response onset.We used a moving window analysis under the assumption that the dynamics would change over time.We also assumed that they would vary relative to different events in the task.
We found significant effects of age for both eigenvalues in the stimulus locked and response locked analyses.Consistent with specificity in age-related changes, we found no significant effects in the cue locked analyses (Fig. 4A,D).The first eigenvalue showed a significant linear effect of age in the stimulus locked condition from 252 to 152 ms before the stimulus, and 15 to 148 (Fig. 4G,I) and 548 to 581 ms after the stimulus (Fig. 4B).The first eigenvalue also showed significant effects in the response locked condition from 248 to 448 ms after the response (Fig. 4C).The second eigenvalue showed significant effects in the stimulus locked analyses from 185 ms before to 248 ms after the stimulus (Fig. 4E,H,J).There were also significant effects of age in the response locked analyses from 219 ms before to 48 ms after the response (Fig. 4H,K), and 215-481 ms after the response (Fig. 4F).Overall, the effects in the second eigenvalue were more widespread in time than the first eigenvalue, but they were locked to specific events.Our model predicts that changes in neural dynamics will be task related.It does not, however, make specific predictions in specific tasks about which epochs would shows changes.However, given the specificity of our effects, our results cannot be accounted for by global changes in the frequency spectrum with age.Given the known difference in the onset of puberty in males and females we performed exploratory analyses including sex as predictor of differences in the eigenvalues.Over all six conditions (2 eigenvalues and three time-locked conditions) and time-points we found a maximum effect of sex for the first eigenvalue in the cue condition with t (305) ¼ 2.59 (p ¼ 0.01), which did not survive correction for multiple comparisons (0.05/6 ¼ 0.0083 Bonferroni corrected.) We next examined correlations between eigenvalues and, behavior including accuracy, reaction times, and the SD of reaction time (Fig. 5).There were no behavioral correlations with the first eigenvalue, across subjects that survived correction for multiple comparisons.In the cue locked analyses, we also found no correlations between eigenvalues and any of the behavioral measures (Fig. 5A), similar to the analyses of age.In the stimulus locked analysis, we found correlations between eigenvalues and reaction times (Fig. 5E) and the SD of the reaction time (Fig. 5H), from 348 to 481 ms after the stimulus (peak at 415 ms, effect estimate 0.233 (SD 0.051) for mean RT and 0.247 (SD 0.056) for RT SD).This is near the average reaction time (Fig. 2C).We also found correlations between the eigenvalues and the reaction time at the time of the response [Fig.5F; from 19 ms before to 48 ms after the response; peak at 15 ms after, effect estimate 0.227 (SD 0.055)] and the SD of the reaction time from 219 ms before until 81 ms and 181-215 ms after the response [Fig.5I, peak at À0.19 ms, effect estimate 0.255 (SD 0.058)].
We also conducted a control analysis, to examine whether possible differences in the timing of the ERPs across age could account for some of our effects.In analyses locked to the response, we controlled for differences in reaction time.However, to further examine the effect of differences in reaction time, we generated a jittered dataset by assuming that the ERP response on single trials was jittered relative to the mean ERP across trials.The jittered dataset was matched for trials and subjects to the actual dataset.For each subject, we generated a new dataset where the ERP in each single trial was equal to the average ERP for the corresponding condition, but with a time offset given by a sample from the SD of reaction times for the corresponding age group (97, 76, and 72 ms for 12, 15, and 18 years old).We then ran these jittered datasets through the same analysis pipeline used to analyze our actual dataset (Fig. 6).Although there were some differences between the groups, they were less pronounced than the differences seen in the actual dataset (Fig. 4).The differences between ages also did not reach statistical significance.Thus, some of the effects seen in the cue or stimulus locked conditions could be driven by variability across trials in processing speed for different subjects, but they likely do not account for a large fraction of our findings.

Discussion
We examined changes in EEG dynamics across adolescence, while subjects performed the Eriksen Flanker task.We fit timevarying linear dynamical systems models to the variability in EEG data, and we examined coefficients (eigenvalues) extracted from the models.Smaller coefficients are consistent with faster decay dynamics, which are consistent with dynamical attractors with steeper basins.We found that the coefficients were inversely correlated with age around the time of stimulus presentation; other features of neural dynamics showed no relations with age, suggesting a level of specificity in developmental effects.We also found positive correlations between reaction time, across subjects, and the model coefficients.Thus, subjects with faster reaction times had smaller coefficients.Overall, the results of these analyses are consistent with the predictions of a computational model, that suggest that adolescent synaptic pruning leads to increased behavioral performance, and changes in neural dynamics (B.B.Averbeck, 2022).
Adolescence is an important developmental period.Performance on basic perceptual tasks achieve adult levels by late childhood (Leat et al., 2009), while cross modal-plasticity between the visual and auditory systems (Cohen et al., 1999) and the ability to learn in some domains (Johnson and Newport, 1989) decreases after late childhood.However, adult performance on tasks requiring executive control does not appear until late adolescence (Luna et al., 2010;Larsen and Luna, 2018).Consistent with these changes in behavior, important changes occur in the brain.Imaging experiments have shown that gray matter volume peaks in late childhood and decrease throughout adolescence, while white matter continues to increase during adolescence (Giedd et al., 1999;Gogtay et al., 2004;Bethlehem et al., 2022).Detailed studies of postmortem tissue have found that the number of excitatory synapses in cortex peaks in late childhood and decreases during adolescence (Huttenlocher, 1979;Huttenlocher and Dabholkar, 1997;Petanjek et al., 2011).Some studies suggest that synaptic pruning in prefrontal cortex occurs later than pruning in unimodal sensory areas (Huttenlocher and Dabholkar, 1997;Elston et al., 2009).Furthermore, myelination increases throughout adolescence (Miller et al., 2012).While it is possible that the changes in gray matter volume measured with structural MR are related to the changes in spines in cortex, spines only make up a small fraction (;1-2%) of the gray matter volume (Rakic et al., 1994), and therefore, this link is not clear.It is possible that synaptic pruning also leads to a decrease in axonal density and dendritic complexity, and these changes contribute to the change in gray matter volume, but this has not, to our knowledge, been investigated.It is also possible that changes in myelination lead to changes in gray matter volume measured with MR, as increased axonal caliber could lead to a change in the MR signal in deep layers of cortex, that is interpreted as a change in gray matter volume (Paus et al., 2008).
Resting state EEG and MEG studies have also found that there are decreases in low-frequency (i.e., u and a) oscillations and increases in high-frequency (i.e., g ) oscillations (Whitford et al., 2007;Segalowitz et al., 2010).These changes have also been found to correlate with changes in gray matter volume (Whitford et al., 2007;Candelaria-Cook et al., 2022).These findings are broadly, although not completely, consistent with our findings.Most importantly, we examined changes in neural dynamics associated with a behavioral task, and not resting state activity.Moreover, we examined specific neural dynamic parameters linked to pruning through our computational model.Of note, our linear dynamical systems model is a first-order (i.e., AR1) low pass filter.In the extreme case, an eigenvalue of 1 in our model would be an integrator, which is a low-pass filter, and an eigenvalue of 0 would be no filtering on a white noise process.Thus, intermediate values correspond to different amounts of low-pass filtering.Also, relations between age and eigenvalues were selective.We found age-related changes in dynamics during specific task epochs, during which specific cognitive operations were unfolding and behaviors were being generated.We found no such correlations during the initial cue period, that unfolds before behaviors manifest.Because there are limited evoked responses during this period, this analysis would be similar to analyzing resting state activity.However, because of the task structure, we did not have sufficient extratask data to estimate resting-state activity.Regardless, our prior work made no predictions linking pruned computational models to resting-state dynamics (B.B.Averbeck, 2022).
We have focused our interpretation, and the corresponding computational modeling, on pruning.Pruning is a well-known computational approach to make large artificial neural networks learn efficiently with less training data (Le Cun et al., 1990;Bishop, 2006).However, it is difficult to measure pruning in children noninvasively, and we examined no direct measure of pruning in the current study.It will be important in future work to attempt to relate changes in synaptic density over adolescence directly to changes in behavior and neural dynamics using animal models, where pruning can be measured directly (Bourgeois et al., 1994;Chung et al., 2017).EEG scalp potentials are thought to be primarily generated by postsynaptic currents form pyramidal neurons, including currents driven by both excitatory and inhibitory synapses.Pruning is thought to affect only excitatory synapses (Petanjek et al., 2011).Other anatomic changes that occur during adolescence, such as myelination, likely shape behavior and corresponding EEG dynamics.Myelination would likely affect conduction speed and perhaps reliability of transmission (Sugio et al., 2023).The effects of myelination on neural dynamics are less clear, as the effects of conduction delays in nonlinear recurrent networks are complex, and their effects on computation are also unclear.Other physiological changes contributing to EEG signal variability, such as changes in skin conductance with development, are also possible confounding factors.In addition, there are numerous other developmental changes.It will, therefore, be important in future work to continue to consider the effects of these alternative mechanisms, and to carry out additional experiments in animal models, that will allow a more detailed examination of the mechanisms giving rise to these changes.A further limitation of our study is the lack of direct measures of brain anatomy changes in the sample.

Limitations
We have conducted our analysis in a principal components (PC) space.PCA extracts components that maximize explained variance, but these components are not necessarily biologically relevant (Delorme et al., 2012).Other approaches to dimensionality reduction could have been used, and future work could compare PCA to other approaches to see which might lead to the most consistent results.
In conclusion, we tested hypotheses from a computational model using changes in task-related variability of EEG signals in a longitudinal sample spanning early to late adolescence.We found, using time-dependent linear dynamical systems models, that the EEG signals showed faster decay dynamics in older subjects.Older subjects also had faster reaction times and increased accuracy, and reaction time correlated with eigenvalues.These results extend computational models that use simulated data.Our data and the associated models suggest that adolescent synaptic pruning leads to increased task performance and steeper attractor basins in recurrent networks trained to perform cognitive tasks.The steeper attractor basins predict the changes in eigenvalues over adolescence in the EEG data.Thus, this suggests that synaptic pruning during adolescence may play a fundamental role in the development of adult competence on tasks requiring executive function, including response inhibition, as measured here with the Eriksen Flanker task.

Figure 1 .
Figure1.Dynamical systems analysis steps.A, Activity from the 104 sensors was first projected into the first two PCs for each subject.This panel shows the time series across a session for a single subject in PC space after z-scoring.B, Distribution of activity for a single subject after z-scoring the first two PCs, showing that data has a zero mean and SD of 1. C, Top, Mean response for a single subject for a single condition in the two PCs.Middle, Single trial evoked response.Bottom, Residual after subtracting the mean response for the corresponding condition from the evoked response.D, Residual activity from a series of trials for a single subject.E, Average and SEM (shaded region) across subjects of the coefficients from the eigenvectors extracted from the A matrix.V is the eigenvector matrix.The first value V 11 corresponds to the coefficient related to the first principal component and the second value V 12 corresponds to the coefficient related to the second principal component.The variable D 1 is the eigenvalue that corresponds to this eigenvector.Bottom, Same analysis applied to the second eigenvector.

Figure 2 .
Figure 2. Task and behavioral results.A, In each trial, subjects were first shown a cue at the center of the screen.After a variable delay, a stimulus was shown, and the subjects were required to indicate whether the central arrow was pointing left or right as quickly and accurately as possible.B, Response accuracy by condition, averaged across age groups.C, Reaction time by condition, averaged across age groups.D, Accuracy by age, averaged across conditions.E, Reaction time by age, averaged across condition.

Figure 3 .
Figure 3. PCA analysis.Data dimensionality was reduced by applying PCA to the grand average evoked responses.A, The first principal component explained 79% of the variance.B, Spatial maps of principal components.C, Subject average ERPs projected into PCA space for incongruent correct and error trials locked to either the stimulus (flan) or the response (resp).

Figure 4 .
Figure 4. Linear effect of age on cue, stimulus, and response locked data.A, Eigenvalue 1 versus time for cue locked data.The fixed effect of the model eigenvalue ; age 1 1|subject is shown as a dashed black line with significant datapoints (p , 0.05/6 cluster corrected) as asterisks.B, Eigenvalue 1 versus time for stimulus locked data.C, Eigenvalue 1 versus time for response locked data.D-F, Same as A-C for eigenvalue 2. G, Bar plots showing the effects of age group on eigenvalue 1 at 15 ms after each task event.H, Same as G for eigenvalue 2. I, Within subject scatter plot of age effect on the first eigenvalue 15 ms after presentation of the stimulus.J. Within subject scatter plot of age effect on second eigenvalue 15 ms after presentation of the stimulus.K, Same as J, 15 ms after the response.

Figure 5 .
Figure 5. Linear effects of the second eigenvalue on subject's accuracy, mean reaction time and SD of reaction time over trials.Red and blue lines show the first and second PC of the ERP averaged across conditions, correspondingly, and dashed line is the time-dependent t statistic.A, Cue locked analysis of correlation between accuracy and eigenvalue.B, Stimulus locked analysis of correlation between accuracy and the second eigenvalue.C, Response locked analysis of correlation between accuracy and eigenvalue.D-F, Same as A-C for correlation between reaction time and the second eigenvalue.G-I, Same as A-C for correlation between the SD of the reaction time and the second eigenvalue.

Figure 6 .
Figure 6.Results of randomizing time offset of evoked response.Panels A-F show eigenvalues (in solid lines, with shaded SE over subjects) and t statistics (in dashed lines) for Cue, Stim, and Response locked analyses of jittered datasets.